! This is a 1D advection example using square initial condition and periodic
! boundary condition for Fromm finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-20: Initial creation.
! - 2024-03-23: Use adv_1d_test_case_mod.

program fromm_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  real, allocatable :: rho(:,:)     ! Tracer density being advected at cell centers
  real, allocatable :: flx(:)       ! Flux at cell interfaces
  logical :: use_rk3 = .false.
  integer, parameter :: ns = 1      ! Stencil width
  integer i
  character(30), parameter :: scheme = 'fromm'

  namelist /params/ nx, nt, dt, use_rk3, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  allocate(flx(1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    ! RK 1st stage
    call fromm(rho(:,old), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - dt / dx * (flx(i) - flx(i-1))
    end do
    call apply_bc(ns, nx, rho(:,new))
    if (use_rk3) then
      ! RK 2nd stage
      call fromm(rho(:,new), flx)
      do i = 1, nx
        rho(i,new) = (3 * rho(i,old) + rho(i,new) - dt / dx * (flx(i) - flx(i-1))) / 4.0d0
      end do
      call apply_bc(ns, nx, rho(:,new))
      ! RK 3st stage
      call fromm(rho(:,new), flx)
      do i = 1, nx
        rho(i,new) = (rho(i,old) + 2 * (rho(i,new) - dt / dx * (flx(i) - flx(i-1)))) / 3.0d0
      end do
      call apply_bc(ns, nx, rho(:,new))
    end if
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  call adv_1d_test_case_final()

  deallocate(rho)
  deallocate(flx)

contains

  subroutine fromm(q, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    real lw, bw
    integer i

    do i = 1, nx
      lw = 0.5d0 * (u * (q(i+1) + q(i)) - dt / dx * u**2 * (q(i+1) - q(i)))
      if (u >= 0) then
        bw = 0.5d0 * (u * (3 * q(i  ) - q(i-1)) - dt / dx * u**2 * (q(i  ) - q(i-1)))
      else
        bw = 0.5d0 * (u * (3 * q(i+1) - q(i+2)) - dt / dx * u**2 * (q(i+2) - q(i+1)))
      end if
      f(i) = 0.5d0 * (lw + bw)
    end do
    call apply_bc(ns, nx, f)

  end subroutine fromm

end program fromm_adv_1d_case
